library(here)
## here() starts at C:/Users/user/Documents/git-stuff/ppm-fire-occurrences
source(here("R", "packages.R"))Non-separable spatio-temporal Poisson point process models for fire occurrences - Companion code
This document replicates the analyses included in the paper entitled Non-separable spatio-temporal Poisson point process models for fire occurrences.
Premilinary steps
Load relevant packages
Download (if necessary) relevant data from the Github Repository. Warning: This code might download a lot of data (around 1.5 GB).
if (!dir.exists(here("data"))) {
dir.create(here("data"))
}
data_names <- c(
"DL_FIRE_J1V-C2_510187.zip", "confini-regioni.zip", "land-use.zip", "NDVI.zip",
"environmental-variables.zip", "INGV-elev.zip"
)
data_paths <- here("data", data_names)
if (!all(file.exists(data_paths))) {
pb_download(
file = data_names,
dest = here("data"),
repo = "agila5/ppm-fire-occurrences",
tag = "v1-data",
overwrite = TRUE,
show_progress = FALSE
)
}
rm(data_names, data_paths)Download (if necessary) cached results from the Github Repository. See the README file for more details.
if (!dir.exists(here("qcache"))) {
dir.create(here("qcache"))
}
qcache_names <- c(
"slope.qs", "NDVI_tidy.qs", "NDVI_tidy_agg_1month.qs", "NDVI_raw.qs", "NDVI_and_landuse.qs", "land_use_tidy.qs", "land_use_tidy_union.qs", "env_var.qs", "elev.qs", "cov_time.qs"
)
qcache_paths <- here("qcache", qcache_names)
if (!all(file.exists(qcache_paths))) {
pb_download(
file = qcache_names,
dest = here("qcache"),
repo = "agila5/ppm-fire-occurrences",
tag = "v1-data",
overwrite = TRUE,
show_progress = FALSE
)
}
rm(qcache_names, qcache_paths)Define several bounding boxes that will be used to create some plots.
define_bb <- function(
xmin,
ymin,
xmax,
ymax,
crs) {
bbox <- st_bbox(
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax),
crs = crs
) |> st_as_sfc()
if (st_crs(bbox) == st_crs(3003)) {
return(bbox)
}
st_transform(bbox, 3003)
}
pantelleria_bbox <- define_bb(
xmin = 11.92586, ymin = 36.73438,
xmax = 12.05684, ymax = 36.83939,
crs = "OGC:CRS84"
)
linosa_bbox <- define_bb(
xmin = 12.84838, ymin = 35.85436,
xmax = 12.88393, ymax = 35.87595,
crs = "OGC:CRS84"
)
lampedusa_bbox <- define_bb(
xmin = 12.51730, ymin = 35.49295,
xmax = 12.63422, ymax = 35.52931,
crs = "OGC:CRS84"
)
palermo_bbox <- define_bb(
xmin = 1871427, xmax = 1890744,
ymin = 4219859, ymax = 4240118,
crs = 3003
)
sicily_mainland_bbox <- define_bb(
xmin = 1807082, ymin = 4041434,
xmax = 2083886, ymax = 4265502,
crs = 3003
)Define also several vectors that will be used to place inset maps into figures.
xrange_pantelleria <- st_bbox(pantelleria_bbox)[c(1, 3)]
xrange_lampedusa <- st_bbox(lampedusa_bbox)[c(1, 3)]
xrange_linosa <- st_bbox(linosa_bbox)[c(1, 3)]
yrange_pantelleria <- st_bbox(pantelleria_bbox)[c(2, 4)]
yrange_lampedusa <- st_bbox(lampedusa_bbox)[c(2, 4)]
yrange_linosa <- st_bbox(linosa_bbox)[c(2, 4)]Section 2
Figure 1
Load relevant data regarding the fires in Italy
fires_italy <- st_read(
paste0("/vsizip/", here("data", "DL_FIRE_J1V-C2_510187.zip")),
quiet = TRUE
) |>
st_transform(3003) |>
mutate(
ACQ_DATETIME = ymd_hm(paste0(ACQ_DATE, " ", ACQ_TIME))
)and borders of the regions
confini_regioni <- st_read(
paste0("/vsizip/", here("data", "confini-regioni.zip"), "/confini-regioni"),
quiet = TRUE
) |>
st_transform(3003)Compute the number of fires in each region
confini_regioni[["counts"]] <- st_intersects(confini_regioni, fires_italy) |> lengths()and plot it
if (!interactive()) {
ggplot(data = confini_regioni) +
geom_sf(aes(fill = counts)) +
scale_fill_continuous_c4a_seq("brewer.oranges")+
theme_minimal() +
theme(
legend.title = element_text(size = 14),
legend.text = element_text(size = 10),
legend.key.height = unit(2, "lines")
) +
labs(fill = "Fire counts")
}Figure 2
Now we need to filter the fires that occurred in Sicily using the bounding box defined by the land-use data object. Therefore, we need to load it
land_use_raw <- st_read(
paste0("/vsizip/", here("data", "land-use.zip"), "/land-use"),
quiet = TRUE
)and apply a series of preprocessing steps to simplify and tidy its structure (following the procedures described in the paper)
land_use_tidy <- qcache(
{
land_use_raw |>
select(Code_18) |>
st_transform(crs = 3003) |> # for spatstat
st_set_agr(c(Code_18 = "constant")) |> # remove warning on "st_cast"
st_cast("POLYGON") |>
st_make_valid() |>
# Merge together areas with the same macro code
mutate(Code_18 = substr(Code_18, 1, 1)) |> # Get the macro code
mutate(
Code_18 = factor(
Code_18,
labels = c(
"Artificial surfaces",
"Agricultural areas",
"Forests",
"Water bodies",
"Water bodies"
)
)
) |>
group_by(Code_18) |>
summarise()
},
name = "land_use_tidy",
cache_dir = here("qcache")
)
## [1] "cached"We need to convert it into tess format for spatstat
owins <-
lapply(
st_geometry(land_use_tidy),
as.owin
) |>
set_names(
land_use_tidy[["Code_18"]]
)
land_use_tess <- tess(tiles = owins); rm(owins)and use its Window attribute to filter the fire points
fires_sicily <- fires_italy[
Window(land_use_tess) |> st_as_sfc() |> st_set_crs(3003),
]There are 10656 events that occurred during 2023 in the region. We can check their monthly temporal distribution and compare with the whole country (Figure 2) as follows:
if (!interactive()) {
bind_rows(
Italy = fires_italy |> st_drop_geometry(),
Sicily = fires_sicily |> st_drop_geometry(),
.id = "ID"
) |>
group_by(ID, month = month(ACQ_DATETIME, label = TRUE)) |>
count() |>
ggplot(aes(x = month, y = n, fill = ID)) +
geom_col(position = position_dodge(), alpha = 0.75) +
geom_text(aes(label = n), fontface = "bold", vjust = 1.5, position = position_dodge(.9), size = 2) +
scale_fill_manual(values = c("orange", "brown")) +
labs(x = "\n Month", y = "Fire Counts\n", fill = "") +
theme_minimal() +
theme(
axis.title = element_text(face = "bold", size = 12)
)
}Figure 3
The following map shows spatial distribution of such events (Figure 3):
land_use_tidy <- land_use_tidy |>
st_set_agr(c(Code_18 = "constant")) |> # remove warning on "st_cast"
st_cast("POLYGON")
land_use_tidy_union <- qcache(
{
land_use_tidy |>
st_geometry() |>
st_union() |>
st_cast("POLYGON")
},
name = "land_use_tidy_union",
cache_dir = here("qcache")
)
## [1] "cached"
idx_eventi_pantelleria <- st_contains(pantelleria_bbox, fires_sicily) |> unlist()
idx_eventi_linosa <- st_contains(linosa_bbox, fires_sicily) |> unlist()
idx_eventi_lampedusa <- st_contains(lampedusa_bbox, fires_sicily) |> unlist()
idx_shape_pantelleria <- st_intersects(pantelleria_bbox, land_use_tidy_union) |> unlist()
shape_pantelleria <- land_use_tidy_union[idx_shape_pantelleria]
idx_shape_linosa <- st_intersects(linosa_bbox, land_use_tidy_union) |> unlist()
shape_linosa <- land_use_tidy_union[idx_shape_linosa]
idx_shape_lampedusa <- st_intersects(lampedusa_bbox, land_use_tidy_union) |> unlist()
shape_lampedusa <- land_use_tidy_union[idx_shape_lampedusa]
shape_everything_else <- land_use_tidy_union[
-c(idx_shape_lampedusa, idx_shape_linosa, idx_shape_pantelleria)
]if (!interactive()) {
# Main plot
main_plot <- ggplot() +
geom_sf(
data = land_use_tidy_union[-c(
idx_shape_lampedusa, idx_shape_linosa, idx_shape_pantelleria
)]
) +
geom_sf(
data = fires_sicily[-c(
idx_eventi_pantelleria, idx_eventi_linosa, idx_eventi_lampedusa
), ]
) +
theme_light() +
theme(panel.background = element_rect(fill = "white"), axis.text = element_text(size = 13.5))
inset_region <- ggplot() +
geom_sf(data = land_use_tidy_union) +
geom_sf(data = lampedusa_bbox |> st_boundary()) +
geom_sf(data = linosa_bbox |> st_boundary()) +
geom_sf(data = pantelleria_bbox |> st_boundary()) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 3),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
inset_pantelleria <- ggplot() +
geom_sf(data = land_use_tidy_union[idx_shape_pantelleria]) +
geom_sf(data = fires_sicily[idx_eventi_pantelleria, ]) +
coord_sf(xlim = c(1760070, 1773730), ylim = c(4068772, 4082023)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 3),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
inset_linosa <- ggplot() +
geom_sf(data = land_use_tidy_union[idx_shape_linosa]) +
geom_sf(data = fires_sicily[idx_eventi_linosa, ]) +
coord_sf(xlim = c(1846520, 1851748), ylim = c(3973711, 3978066)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 3),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
inset_lampedusa <- ggplot() +
geom_sf(data = land_use_tidy_union[idx_shape_lampedusa]) +
geom_sf(data = fires_sicily[idx_eventi_lampedusa, ]) +
coord_sf(xlim = c(1818045, 1830515), ylim = c(3932694, 3938490)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 3),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
ggplot() +
coord_equal(xlim = c(0, 24), ylim = c(0, 19), expand = FALSE) +
annotation_custom(
ggplotGrob(main_plot),
xmin = 0, xmax = 24, ymin = 0, ymax = 19
) +
annotation_custom(
ggplotGrob(inset_region),
xmin = 1.5, xmax = 7, ymin = 12.5, ymax = 18
) +
annotation_custom(
ggplotGrob(inset_pantelleria),
xmin = 2, xmax = 5, ymin = 2, ymax = 5
) +
annotation_custom(
ggplotGrob(inset_lampedusa),
xmin = 5.75,
xmax = 5.75 + 3 * diff(xrange_lampedusa) / diff(xrange_pantelleria),
ymin = 2.6,
ymax = 2.6 + 5 * diff(yrange_lampedusa) / diff(yrange_pantelleria)
) +
annotation_custom(
ggplotGrob(inset_linosa),
xmin = 9.1,
xmax = 9.1 + 3 * diff(xrange_linosa) / diff(xrange_pantelleria),
ymin = 3,
ymax = 3 + 5 * diff(yrange_linosa) / diff(yrange_pantelleria)
) +
annotate(
"segment",
x = 3.5, y = 5.1,
xend = 2.825, yend = 10.2,
lineend = "round",
linewidth = 1
) +
annotate(
"segment",
x = 2.7, y = 11,
xend = 2.25, yend = 14.3,
arrow = arrow(),
lineend = "round",
linewidth = 1
) +
annotate(
"segment",
x = 7.1, y = 4.3,
xend = 4.355, yend = 10,
lineend = "round",
linewidth = 1
) +
annotate(
"segment",
x = 4.05, y = 10.7,
xend = 3.1, yend = 12.65,
arrow = arrow(),
lineend = "round",
linewidth = 1,
) +
annotate(
"segment",
x = 9.5, y = 4,
xend = 7, yend = 7.75,
lineend = "round",
linewidth = 1,
) +
annotate(
"segment",
x = 4.85, y = 11,
xend = 3.4, yend = 13.2,
arrow = arrow(),
lineend = "round",
linewidth = 1,
) +
labs(x = "", y = "") +
theme(
panel.background = element_rect(fill = "white"),
axis.text = element_blank(),
axis.ticks = element_blank()
)
}Figure 4
Spatio-temporal monthly distribution of the events occurred in the mainland. First we need to define such mainland
idx_mainland <- which.max(st_area(land_use_tidy_union))
mainland <- land_use_tidy_union[idx_mainland]; rm(idx_mainland)and then we can plot such events
if (!interactive()) {
fires_sicily |>
st_filter(mainland) |>
mutate(month = month(ACQ_DATETIME, label = TRUE, abbr = FALSE)) |>
ggplot() +
geom_sf(data = mainland) +
geom_sf() +
facet_wrap(~month) +
theme_light() +
scale_x_continuous(breaks = c(12.5, 13.5, 14.5, 15.5)) +
theme(
strip.text = element_text(face = "bold", colour = "black")
)
}Figure 5
We need to define a series of indices that will be used in the inset maps
idx_pantelleria <- st_intersects(pantelleria_bbox, land_use_tidy) |> unlist()
idx_linosa <- st_intersects(linosa_bbox, land_use_tidy) |> unlist()
idx_lampedusa <- st_intersects(lampedusa_bbox, land_use_tidy) |> unlist()
land_use_palermo <- st_intersection(land_use_tidy, palermo_bbox)and then we can visualise the land use variable
if (!interactive()) {
mainland_plot <- ggplot() +
geom_sf(
data = land_use_tidy[-c(
idx_lampedusa, idx_linosa, idx_pantelleria
), ],
aes(fill = Code_18)
) +
geom_sf(
data = st_boundary(palermo_bbox),
linewidth = 1
) +
theme_light() +
scale_fill_manual(
values = c(
"Artificial surfaces" = "#a50000",
"Agricultural areas" = "#e49703",
"Forests" = "#287201",
"Water bodies" = "#bdeafe"
)
) +
theme(
axis.text = element_text(size = 13.5),
legend.text = element_text(size = 15),
plot.title = element_text(face = "bold", size = 18, hjust = 0.5)
) +
labs(fill = "", title = "Land usage")
palermo_zoom <- ggplot() +
geom_sf(
data = land_use_palermo,
aes(fill = Code_18),
show.legend = FALSE
) +
geom_sf(
data = st_boundary(palermo_bbox),
linewidth = 1.5
) +
theme_light() +
coord_sf(datum = st_crs(3003)) +
scale_fill_manual(
values = c(
"Artificial surfaces" = "#a50000",
"Agricultural areas" = "#e49703",
"Forests" = "#287201",
"Water bodies" = "#bdeafe"
)
) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", colour = NA),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
pantelleria_zoom <- ggplot() +
geom_sf(
data = land_use_tidy[idx_pantelleria, ] |> st_transform(4326),
aes(fill = Code_18),
show.legend = FALSE
) +
geom_sf(
data = st_boundary(land_use_tidy_union[pantelleria_bbox, ]),
linewidth = 1
) +
theme_light() +
scale_fill_manual(
values = c(
"Artificial surfaces" = "#a50000",
"Agricultural areas" = "#e49703",
"Forests" = "#287201",
"Water bodies" = "#bdeafe"
)
) +
theme(
panel.border = element_rect(fill = NA, linewidth = 2, colour = "black"),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
panel.grid = element_blank()
)
lampedusa_zoom <- ggplot() +
geom_sf(
data = land_use_tidy[idx_lampedusa, ] |> st_transform(4326),
aes(fill = Code_18),
show.legend = FALSE
) +
geom_sf(
data = st_boundary(land_use_tidy_union[lampedusa_bbox, ]),
linewidth = 1
) +
theme_light() +
scale_fill_manual(
values = c(
"Artificial surfaces" = "#a50000",
"Agricultural areas" = "#e49703",
"Forests" = "#287201",
"Water bodies" = "#bdeafe"
)
) +
theme(
panel.border = element_rect(fill = NA, linewidth = 2, colour = "black"),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
panel.grid = element_blank()
)
linosa_zoom <- ggplot() +
geom_sf(
data = land_use_tidy[idx_linosa, ] |> st_transform(4326),
aes(fill = Code_18),
show.legend = FALSE
) +
geom_sf(
data = st_boundary(land_use_tidy_union[linosa_bbox, ]),
linewidth = 1
) +
theme_light() +
scale_fill_manual(
values = c(
"Artificial surfaces" = "#a50000",
"Agricultural areas" = "#e49703",
"Forests" = "#287201",
"Water bodies" = "#bdeafe"
)
) +
theme(
panel.border = element_rect(fill = NA, linewidth = 1, colour = "black"),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
panel.grid = element_blank()
)
ggplot() +
coord_equal(xlim = c(0, 24), ylim = c(0, 20), expand = FALSE) +
annotation_custom(
ggplotGrob(mainland_plot),
xmin = 0, xmax = 24, ymin = 0, ymax = 20
) +
annotation_custom(
ggplotGrob(palermo_zoom),
xmin = 1.8, xmax = 6, ymin = 12.15, ymax = 16.65
) +
annotation_custom(
ggplotGrob(pantelleria_zoom),
xmin = 2, xmax = 5, ymin = 4, ymax = 7
) +
annotation_custom(
ggplotGrob(lampedusa_zoom),
xmin = 5.25, xmax = 5.25 + 3 * diff(xrange_lampedusa) / diff(xrange_pantelleria),
ymin = 4.6, ymax = 4.6 + 5 * diff(yrange_lampedusa) / diff(yrange_pantelleria)
) +
annotation_custom(
ggplotGrob(linosa_zoom),
xmin = 8.15, xmax = 8.15 + 3 * diff(xrange_linosa) / diff(xrange_pantelleria),
ymin = 5, ymax = 5 + 5 * diff(yrange_linosa) / diff(yrange_pantelleria)
) +
annotate(
"segment",
x = 6, y = 14.75,
xend = 8, yend = 12.75,
arrow = arrow(),
lineend = "round",
linewidth = 1,
) +
theme(
panel.background = element_rect(fill = "white"),
axis.text = element_blank(),
axis.ticks = element_blank()
) +
labs(x = "", y = "")
}Figure 6
Graphical representation of the moving grid for Horn’s Algorithm
if (!interactive()) {
horn_grid <- st_make_grid(
cellsize = c(1, 1), offset = c(0, 0), n = c(3, 3)
)
ggplot() +
geom_sf(
data = horn_grid[5],
fill = grey(0.85)
) +
geom_sf(data = st_boundary(horn_grid), linewidth = 1) +
geom_sf_text(
data = st_centroid(horn_grid),
label = c( # NB: The grid is specified in reverse order
"Alt[3]", "Alt[4]", "Alt[5]",
"Alt[2]", "", "Alt[6]",
"Alt[1]", "Alt[8]", "Alt[7]"
),
parse = TRUE,
size = 13,
fontface = "bold"
) +
labs(x = "Longitude", y = "Latitude") +
theme(
panel.background = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_text(size = 40, face = "bold")
)
}Figure 7a
Now we need to focus on the Elevation in Sicily. First, we define an auxiliary function
generate_orig_tif <- function(paths, tmp_dir) {
unzip(paths, exdir = tmp_dir)
tif_paths <- list.files(
tmp_dir, pattern = "\\.tif", recursive = TRUE,
full.names = TRUE
)
tifs <- lapply(tif_paths, read_stars)
out <- do.call(st_mosaic, tifs); invisible(gc())
out
}and then read-in the data
elev <- qcache(
{
tmp_dir <- tempdir()
orig_tif <- generate_orig_tif(here("data", "INGV-elev.zip"), tmp_dir)
# Clean files
unlink(paste0(tmp_dir, "INGV-elev"), recursive = TRUE)
invisible(gc())
tif <- orig_tif |>
st_downsample(n = c(10, 10)) |>
st_warp(crs = 3003)
# Clean files
unlink(list.files(tmp_dir, pattern = "\\.tif", full.names = TRUE))
rm(orig_tif)
invisible(gc())
tif <- tif[st_bbox(land_use_tidy_union)]
tif[land_use_tidy_union]
},
name = "elev",
cache_dir = here("qcache")
)
## [1] "cached"
invisible(gc())Now we can represent it
elev_mainland <- (elev |> st_downsample(c(5, 5)))[shape_everything_else]
elev_pantelleria <- elev[shape_pantelleria]
elev_lampedusa <- elev[shape_lampedusa]
elev_linosa <- elev[shape_linosa]if (!interactive()) {
elev_mainplot <- ggplot() +
geom_sf(
aes(fill = e41005_s10.tif),
data = elev_mainland |> st_as_sf(),
linewidth = NA
) +
theme_light() +
scale_fill_gradientn(
breaks = c(0, 500, 1000, 2000, 3000),
colours = terrain.colors(100),
trans = modulus_trans(p = 0.8)
) +
theme(
axis.text = element_text(size = 19.5),
legend.text = element_text(size = 19.5),
legend.key.height = unit(3, "lines"),
plot.title = element_text(face = "bold", size = 26, hjust = 0.5)
) +
labs(fill = "", title = "Altitude (m)")
elev_pantelleria_plot <- ggplot() +
geom_stars(data = elev_pantelleria |> st_warp(crs = 4326), na.action = na.omit, show.legend = FALSE) +
coord_sf(crs = 4326) +
scale_fill_gradientn(colours = terrain.colors(100), trans = modulus_trans(p = 0.8), limits = c(0, 3000)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 3),
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
elev_lampedusa_plot <- ggplot() +
geom_stars(
data = elev_lampedusa |>
st_warp(crs = 4326),
na.action = na.omit,
show.legend = FALSE
) +
coord_sf(crs = 4326) +
scale_fill_gradientn(colours = terrain.colors(100), trans = modulus_trans(p = 0.7), limits = c(0, 3000)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 3),
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
elev_linosa_plot <- ggplot() +
geom_stars(
data = elev_linosa |>
st_warp(crs = 4326),
na.action = na.omit,
show.legend = FALSE
) +
coord_sf(crs = 4326) +
scale_fill_gradientn(colours = terrain.colors(100), trans = modulus_trans(p = 0.7), limits = c(0, 3000)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 2),
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
ggplot() +
coord_equal(xlim = c(0, 24), ylim = c(0, 20), expand = FALSE) +
annotation_custom(
ggplotGrob(elev_mainplot),
xmin = 0, xmax = 24, ymin = 0, ymax = 20
) +
annotation_custom(
ggplotGrob(elev_pantelleria_plot),
xmin = 2.5, xmax = 5.5, ymin = 3.15, ymax = 6.15
) +
annotation_custom(
ggplotGrob(elev_lampedusa_plot),
xmin = 6.25, xmax = 6.25 + 3 * diff(xrange_lampedusa) / diff(xrange_pantelleria),
ymin = 3.8, ymax = 3.8 + 3 * diff(yrange_lampedusa) / diff(yrange_pantelleria)
) +
annotation_custom(
ggplotGrob(elev_linosa_plot),
xmin = 9.5, xmax = 9.5 + 3 * diff(xrange_linosa) / diff(xrange_pantelleria),
ymin = 4.2, ymax = 4.2 + 3 * diff(yrange_linosa) / diff(yrange_pantelleria)
) +
theme(panel.background = element_rect(fill = "white"))
}Figure 7b
As we mentioned in the paper, the slope is computed using the GDAL DEM tools
slope <- qcache(
{
tmp_dir <- tempdir()
orig_tif <- generate_orig_tif(here("data", "INGV-elev.zip"), tmp_dir)
unlink(paste0(tmp_dir, "INGV-elev"), recursive = TRUE)
invisible(gc())
temp_tif <- tempfile(fileext = ".tif")
write_stars(orig_tif, temp_tif, progress = FALSE)
temp_slope <- tempfile(fileext = ".tif")
gdaldem("slope", temp_tif, temp_slope)
slope <- read_stars(temp_slope)
slope <- slope |> st_downsample(c(10, 10)) |> st_warp(crs = 3003)
slope <- slope[st_bbox(land_use_tidy_union)]
slope <- slope[land_use_tidy_union]
unlink(list.files(tmp_dir, pattern = "\\.tif", full.names = TRUE))
rm(orig_tif, temp_tif, temp_slope)
invisible(gc())
slope
},
name = "slope",
cache_dir = here("qcache")
)
## [1] "cached"As before, we can subset the stars object
slope_mainland <- (slope |> st_downsample(c(5, 5)))[shape_everything_else]
slope_pantelleria <- slope[shape_pantelleria]
slope_lampedusa <- slope[shape_lampedusa]
slope_linosa <- slope[shape_linosa]and now we can plot it
if (!interactive()) {
slope_mainplot <- ggplot() +
geom_sf(
aes(fill = value),
data = slope_mainland |> st_as_sf() |> setNames(c("value", "geometry")),
linewidth = NA
) +
theme_light() +
scale_fill_gradientn(
colours = terrain.colors(100),
trans = modulus_trans(p = 0.8)
) +
theme(
axis.text = element_text(size = 19.5),
legend.text = element_text(size = 19.5),
legend.key.height = unit(3, "lines"),
plot.title = element_text(face = "bold", size = 26, hjust = 0.5)
) +
labs(fill = "", title = "Slope (°)")
slope_pantelleria_plot <- ggplot() +
geom_stars(data = slope_pantelleria |> st_warp(crs = 4326), na.action = na.omit, show.legend = FALSE) +
coord_sf(crs = 4326) +
scale_fill_gradientn(colours = terrain.colors(100), trans = modulus_trans(p = 0.8), limits = c(0, 80)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 3),
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
slope_lampedusa_plot <- ggplot() +
geom_stars(data = slope_lampedusa |> st_warp(crs = 4326), na.action = na.omit, show.legend = FALSE) +
coord_sf(crs = 4326) +
scale_fill_gradientn(colours = terrain.colors(100), trans = scales::modulus_trans(p = 0.7), limits = c(0, 80)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 3),
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
slope_linosa_plot <- ggplot() +
geom_stars(data = slope_linosa |> st_warp(crs = 4326), na.action = na.omit, show.legend = FALSE) +
coord_sf(crs = 4326) +
scale_fill_gradientn(colours = terrain.colors(100), trans = scales::modulus_trans(p = 0.7), limits = c(0, 80)) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", linewidth = 2),
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
plot.margin = margin(0, 0, 0, 0),
panel.grid = element_blank()
)
ggplot() +
coord_equal(xlim = c(0, 24), ylim = c(0, 20), expand = FALSE) +
annotation_custom(
ggplotGrob(slope_mainplot),
xmin = 0, xmax = 24, ymin = 0, ymax = 20
) +
annotation_custom(
ggplotGrob(slope_pantelleria_plot),
xmin = 2.5, xmax = 5.5, ymin = 3.15, ymax = 6.15
) +
annotation_custom(
ggplotGrob(slope_lampedusa_plot),
xmin = 6.25, xmax = 6.25 + 3 * diff(xrange_lampedusa) / diff(xrange_pantelleria),
ymin = 3.8, ymax = 3.8 + 3 * diff(yrange_lampedusa) / diff(yrange_pantelleria)
) +
annotation_custom(
ggplotGrob(slope_linosa_plot),
xmin = 9.5, xmax = 9.5 + 3 * diff(xrange_linosa) / diff(xrange_pantelleria),
ymin = 4.2, ymax = 4.2 + 3 * diff(yrange_linosa) / diff(yrange_pantelleria)
) +
theme(panel.background = element_rect(fill = "white"))
}Figure 8
The NDVI data can be loaded as follows
NDVI_raw <- qcache(
{
tmp_dir <- tempdir()
unzip(here("data", "NDVI.zip"), exdir = tmp_dir)
files <- list.files(
paste0(tmp_dir, "/NDVI"),
pattern = "\\.nc$",
full.names = TRUE
)
# Currently, the structure returned by read_stars is an array with
# dimension 1253 x 1116 and 36 attributes (1 for each NDVI in a 10 days
# interval). We need to convert it to a 1253 x 1116 x 36 array with 1
# attribute. See below (NDVI_tidy).
out <- read_stars(files, quiet = TRUE, proxy = FALSE)
unlink(paste0(tmp_dir, "/NDVI"), recursive = TRUE)
out
},
name = "NDVI_raw",
cache_dir = here("qcache")
)
## [1] "cached"As mentioned in the code chunk, the NDVI_raw object needs a bit of preprocessing to convert it into a usable format
# The online docs say that the NDVI measurements are shared approximately
# every 10 days
day_jan <- c(
as.Date("2023-01-01"), as.Date("2023-01-11"), as.Date("2023-01-21")
)
# The following code is taken from the first vignette of stars package.
NDVI_tidy <- qcache(
{
merge(NDVI_raw) |>
setNames("NDVI") |>
st_set_dimensions(
3,
values = c(
day_jan,
day_jan + months(1),
day_jan + months(2),
day_jan + months(3),
day_jan + months(4),
day_jan + months(5),
day_jan + months(6),
day_jan + months(7),
day_jan + months(8),
day_jan + months(9),
day_jan + months(10),
day_jan + months(11)
),
names = "time"
) |>
st_warp(crs = 3003, cellsize = c(500, 500))
},
name = "NDVI_tidy",
cache_dir = here("qcache")
)
## [1] "cached"
rm(day_jan)
invisible(gc())We can temporally aggregate the NDVI values at the monthly level
NDVI_tidy_agg_1month <- qcache(
{
aggregate(NDVI_tidy, by = "months", FUN = mean)[mainland]
},
name = "NDVI_tidy_agg_1month",
cache_dir = here("qcache")
)
## [1] "cached"and show them
if (!interactive()) {
ggplot() +
geom_stars(data = NDVI_tidy_agg_1month, downsample = c(1, 1, 0)) +
facet_wrap(~time) +
scale_fill_gradient2(
low = "#8d5e0b",
mid = "#e6e91a",
high = "#045b00",
midpoint = 0.5,
limits = c(-0.08, 0.92),
na.value = "lightblue"
) +
theme_light() +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
panel.grid = element_blank(),
panel.background = element_rect(fill = "lightblue"),
strip.text = element_text(face = "bold", size = 13),
legend.title = element_text(size = 13),
legend.text = element_text(size = 10)
)
}Figure 9
The following code computes the average NDVI for each land use type
NDVI_and_landuse <- qcache(
{
st_warp(
NDVI_tidy,
st_as_stars(st_bbox(NDVI_tidy), dx = 1500, dy = 1500)
) |>
aggregate(
by = land_use_tidy |> group_by(Code_18) |> summarise(),
FUN = mean,
na.rm = TRUE
)
},
name = "NDVI_and_landuse",
cache_dir = here("qcache")
)
## [1] "cached"and we can plot it as follows
if (!interactive()) {
NDVI_and_landuse |>
st_as_sf() |>
mutate(idx = c("Artificial surfaces", "Agricultural areas", "Forests", "Water bodies")) |>
st_drop_geometry() |>
pivot_longer(!idx, names_to = "date", names_transform = list(
date = as.Date
)) |>
ggplot() +
geom_line(aes(x = date, y = value, col = idx, linetype = idx), linewidth = 1.3) +
theme_light() +
labs(x = "", y = "Average NDVI", linetype = "", col = "")
}Figure 10
First, we need to define a few auxiliary functions:
normalize <- function(x, max = 1) {
x <- as.numeric(x)
(x - min(x)) / (max(x) - min(x)) * max
}
# The following function is used to create the space-time representation
# of the environmental variables. See below.
tf = function(x, y, w = .75, h = .33) {
x2 = x * w + y * (1 - w)
y2 = y * h
x2[length(x2) + 1] = x2[1]
y2[length(y2) + 1] = y2[1]
list(x = x2, y = y2)
}Next we load the environmental variables downloaded from ERA5
env_var <- qcache(
{
tmp_dir <- tempdir()
unzip(here("data", "environmental-variables.zip"), exdir = tmp_dir)
nc_files <- list.files(
paste0(tmp_dir, "/environmental-variables"),
pattern = "\\.nc$",
full.names = TRUE
)
img_combined <- lapply(
nc_files,
function(x, var) {
img <- read_ncdf(x, var = var)
if (grepl("october.nc", x, fixed = TRUE)) {
img <- img[, , , 2, , drop = TRUE]
}
img
},
var = c("u10", "v10", "d2m", "t2m", "skt", "stl1", "stl2", "stl3", "stl4", "sp", "tp")
) |>
do.call(c, args = _)
unlink(nc_files); rm(nc_files); invisible(gc())
st_crs(img_combined) <- "OGC:CRS84"
st_warp(img_combined, crs = 3003)
},
name = "env_var",
cache_dir = here("qcache")
)
## [1] "cached"
env_var <- env_var[mainland]; invisible(gc()) # Subset mainland
env_var <- aggregate(env_var, max, by = "1 day"); invisible(gc()) # Take daily maximumNow we can replicate the space-time representation of the environmental variables. We start from Surface Pressure
# Subset the surface pressure data and build the palette
sp_2months <- aggregate(
env_var["sp", ],
mean,
by = "2 months",
na.rm = TRUE
)
sp_palette <- col_numeric(
viridis::turbo(15),
domain = sp_2months$sp |> as.vector() |> range(na.rm = TRUE)
)if (!interactive()) {
grid.newpage()
vp1 <- viewport(
x = unit(0.35, "npc"),
y = unit(0.35, "npc"),
width = unit(0.8, "npc"),
height = unit(0.8, "npc")
)
vp2 <- viewport(
x = unit(0.9, "npc"),
y = unit(0.35, "npc"),
width = unit(0.25, "npc"),
height = unit(0.6, "npc")
)
pushViewport(vp1)
w <- 0.74
h <- 0.175
for (time_period in 1:6) {
(sp_2months[, time_period, drop = TRUE]) |>
st_as_sf() -> sp_2months_sf
poly_cols <- sp_palette(sp_2months_sf$sp)
poly_coords <- sp_2months_sf |>
st_geometry() |>
st_coordinates()
poly_coords[, 1] <- normalize(poly_coords[, 1])
poly_coords[, 2] <- normalize(poly_coords[, 2])
xy_trans <- tf(poly_coords[, 1], poly_coords[, 2], w = w, h = h)
for (i in unique(poly_coords[, 4])) {
idx <- which(poly_coords[, 4] == i)
grid.polygon(
x = xy_trans$x[idx] + 0.075,
y = xy_trans$y[idx] + 0.1 + (time_period - 1) * 0.9 / 12,
gp = gpar(fill = poly_cols[i], col = NA, alpha = 1)
)
}
sicily <- sp_2months_sf |>
st_union() |>
st_geometry()
sicily_coords <- sicily |> st_coordinates()
sicily_coords[, 1] <- normalize(sicily_coords[, 1])
sicily_coords[, 2] <- normalize(sicily_coords[, 2])
xy_trans <- tf(sicily_coords[, 1], sicily_coords[, 2], w = w, h = h)
grid.polygon(
x = xy_trans$x + 0.075,
y = xy_trans$y + 0.1 + (time_period - 1) * 0.9 / 12,
gp = gpar(fill = NA, col = "black", lwd = 2)
)
}
grid.text(
"Jan-Feb",
x = unit(0.15, "npc"),
y = unit(0.19, "npc"),
gp = gpar(fontsize = 13, fontface = "bold")
)
grid.text(
"Nov-Dec",
x = unit(0.15, "npc"),
y = unit(0.59, "npc"),
gp = gpar(fontsize = 13, fontface = "bold")
)
grid.text(
"Surface Pressure (Pa)",
x = unit(0.65, "npc"),
y = unit(0.7, "npc"),
gp = gpar(fontsize = 20, fontface = "bold")
)
grid.segments(
x0 = unit(0.15, "npc"),
x1 = unit(0.15, "npc"),
y0 = unit(0.21, "npc"),
y1 = unit(0.565, "npc"),
gp = gpar(lwd = 2),
arrow = arrow()
)
popViewport()
pushViewport(vp2)
labels <- seq(88099, 102800, length.out = 5) |> pretty()
grid.points(
x = unit(rep(0.3, 5), "npc"),
y = unit(seq(0.2, 0.6, length.out = 5), "npc"),
pch = 21,
size = unit(1.5, "char"),
gp = gpar(fill = sp_palette(seq(89000, 102000, length.out = 5)))
)
grid.text(
label = labels,
x = unit(rep(0.55, 5), "npc"),
y = unit(seq(0.2, 0.6, length.out = 5), "npc"),
gp = gpar(fontsize = 15, fontface = "bold")
)
}and then focus on Skin Temperature
stl2_2months <- aggregate(
env_var["stl2", ],
mean,
by = "2 months",
na.rm = TRUE
)
stl2_2months$stl2 <- stl2_2months$stl2 - 273.15
stl2_palette <- col_numeric(
viridis::turbo(15),
domain = stl2_2months$stl2 |> as.vector() |> range(na.rm = TRUE)
)if (!interactive()) {
grid.newpage()
vp1 <- viewport(
x = unit(0.35, "npc"),
y = unit(0.35, "npc"),
width = unit(0.8, "npc"),
height = unit(0.8, "npc")
)
vp2 <- viewport(
x = unit(0.9, "npc"),
y = unit(0.35, "npc"),
width = unit(0.25, "npc"),
height = unit(0.6, "npc")
)
pushViewport(vp1)
w <- 0.74
h <- 0.175
for (time_period in 1:6) {
(stl2_2months[, time_period, drop = TRUE]) |>
st_as_sf() -> stl2_2months_sf
poly_cols <- stl2_palette(stl2_2months_sf$stl2)
poly_coords <- stl2_2months_sf |>
st_geometry() |>
st_coordinates()
poly_coords[, 1] <- normalize(poly_coords[, 1])
poly_coords[, 2] <- normalize(poly_coords[, 2])
xy_trans <- tf(poly_coords[, 1], poly_coords[, 2], w = w, h = h)
for (i in unique(poly_coords[, 4])) {
idx <- which(poly_coords[, 4] == i)
grid.polygon(
x = xy_trans$x[idx] + 0.075,
y = xy_trans$y[idx] + 0.1 + (time_period - 1) * 0.9 / 12,
gp = gpar(fill = poly_cols[i], col = NA, alpha = 1)
)
}
sicily <- stl2_2months_sf |>
st_union() |>
st_geometry()
sicily_coords <- sicily |> st_coordinates()
sicily_coords[, 1] <- normalize(sicily_coords[, 1])
sicily_coords[, 2] <- normalize(sicily_coords[, 2])
xy_trans <- tf(sicily_coords[, 1], sicily_coords[, 2], w = w, h = h)
grid.polygon(
x = xy_trans$x + 0.075,
y = xy_trans$y + 0.1 + (time_period - 1) * 0.9 / 12,
gp = gpar(fill = NA, col = "black", lwd = 2)
)
}
grid.text(
"Jan-Feb",
x = unit(0.15, "npc"),
y = unit(0.19, "npc"),
gp = gpar(fontsize = 13, fontface = "bold")
)
grid.text(
"Nov-Dec",
x = unit(0.15, "npc"),
y = unit(0.59, "npc"),
gp = gpar(fontsize = 13, fontface = "bold")
)
grid.text(
"Skin Temperature (Celsius)",
x = unit(0.65, "npc"),
y = unit(0.7, "npc"),
gp = gpar(fontsize = 20, fontface = "bold")
)
grid.segments(
x0 = unit(0.15, "npc"),
x1 = unit(0.15, "npc"),
y0 = unit(0.21, "npc"),
y1 = unit(0.565, "npc"),
gp = gpar(lwd = 2),
arrow = arrow()
)
popViewport()
pushViewport(vp2)
labels <- seq(5, 35, length.out = 5) |> pretty()
grid.points(
x = unit(rep(0.3, 7), "npc"),
y = unit(seq(0.1, 0.6, length.out = 7), "npc"),
pch = 21,
size = unit(1.5, "char"),
gp = gpar(fill = stl2_palette(seq(6.5, 33, length.out = 7)))
)
grid.text(
label = labels,
x = unit(rep(0.5, 7), "npc"),
y = unit(seq(0.1, 0.6, length.out = 7), "npc"),
gp = gpar(fontsize = 18, fontface = "bold")
)
}Section 4
Figure 11
Now we replicate the plots regarding the temporal distribution of some environmental variables. First we need to reload the data (since we modified it to create the previous plot) and adjust its format slightly:
cov_time <- qcache(
{
tmp_dir <- tempdir()
unzip(here("data", "environmental-variables.zip"), exdir = tmp_dir)
nc_files <- list.files(
paste0(tmp_dir, "/environmental-variables"),
pattern = "\\.nc$",
full.names = TRUE
)
nome_file <- c(
"january.nc", "february.nc", "march.nc", "april.nc",
"may.nc", "june.nc", "july.nc",
"august.nc", "september.nc", "october.nc", "november.nc",
"december.nc"
)
lista <- list()
for (i in nc_files) {
k <- which(nome_file == basename(i))
suppressMessages({
prova <- read_ncdf(i)
})
st_crs(prova) <- 4326
prova <- st_warp(prova, crs = 3003)
prova <- prova[mainland]
if (basename(i) == "october.nc") {
prova2 <- aggregate(prova[, , , 2, , drop = TRUE], max, by = "1 day")[, -1] |>
st_apply(c("time"), function(z) max(z, na.rm = TRUE))
} else {
prova2 <- aggregate(prova[, , , ], max, by = "1 day") |>
st_apply(c("time"), function(z) max(z, na.rm = TRUE))
}
prova2 <- as.data.frame(prova2)
prova2$month <- rep(k, nrow(prova2))
prova2$day <- 1:nrow(prova2)
lista[[k]] <- prova2
rm(prova, prova2)
invisible(gc())
}
rm(nome_file, i, k)
unlink(nc_files); rm(nc_files); invisible(gc())
do.call(rbind, lista)
},
name = "cov_time",
cache_dir = here("qcache")
)
## [1] "cached"if (exists("lista")) {
rm(lista)
}We need to convert Kelvin to Celsius
cov_time$d2m <- cov_time$d2m - 273.15
cov_time$t2m <- cov_time$t2m - 273.15
cov_time$skt <- cov_time$skt - 273.15
cov_time$stl1 <- cov_time$stl1 - 273.15
cov_time$stl2 <- cov_time$stl2 - 273.15
cov_time$stl3 <- cov_time$stl3 - 273.15
cov_time$stl4 <- cov_time$stl4 - 273.15and now we are ready to plot
if (!interactive()) {
par(mfrow = c(2, 2))
plot(
cov_time$time, cov_time$u10, type = "l", col = 5,
xlab = "", ylab = "m/s", main = "Wind Speed", lwd = 2
)
lines(cov_time$time, cov_time$v10, col = 6, lwd = 2)
legend("top", legend = c("u10", "v10"), col = c(5, 6), lty = 1, cex = 0.8, lwd = 2)
plot(
cov_time$time, cov_time$d2m, type = "l", col = 4,
xlab = "", ylab = "Celsius", main = "Temperatures", lwd = 2,
ylim = c(2, 60)
)
lines(cov_time$time, cov_time$t2m, col = 5, lwd = 2)
lines(cov_time$time, cov_time$skt, col = 6, lwd = 2)
lines(cov_time$time, cov_time$stl1, col = 7, lwd = 2)
lines(cov_time$time, cov_time$stl2, col = 8, lwd = 2)
lines(cov_time$time, cov_time$stl3, col = 9, lwd = 2)
lines(cov_time$time, cov_time$stl4, col = 10, lwd = 2)
legend("topleft", legend=c("d2m", "t2m", "skt", "stl1", "stl2", "stl3", "stl4"),
col=c(4:10), lty=1, cex=0.8, lwd = 2)
plot(
cov_time$time, cov_time$tp, type = "l", col = 4,
xlab = "", ylab = "m", main = "Precipitations", lwd = 2
)
plot(
cov_time$time, cov_time$sp, type = "l", col = 4,
xlab = "", ylab = "Pa", main = "Surface Pressure", lwd = 2
)
}Figure 12
Next we focus on the pairs plot
panel.hist <- function(x, ...)
{
usr <- par("usr")
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col = "grey", ...)
}
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
if (!interactive()) {
pairs(
cov_time[,c(2:12)], upper.panel = panel.cor, diag.panel = panel.hist,
lower.panel = panel.smooth, pch = "."
)
}For the estimation of the statistical model detailed in Section 3 of the manuscript, we need to define a spatio-temporal structure that contains the points
stp_time <- normalize(fires_sicily[["ACQ_DATETIME"]], 364)
stp0 <- stp(cbind(st_coordinates(fires_sicily) / 1000, stp_time))and load a set of randomly defined (spatio-temporal) dummy points
dum1 <- qread(here("qcache", "dum1.qs"))Then we need to extract or compute the value of each covariate for these observed or dummy points. For the purely spatial covariates we can simply used a lookup table as defined in spatstat:
elev_im <- spatstat.geom::rescale(as.im(elev), 1000); rm(elev); invisible(gc())
slope_im <- spatstat.geom::rescale(as.im(slope), 1000); rm(slope); invisible(gc())
land_use_f <- as.function(land_use_tess)
covs0_spatial <- qcache(
{
rbind(
cbind(
land_use_f(stp0$df$x * 1000, stp0$df$y * 1000),
slope_im[list(x = stp0$df$x, y = stp0$df$y), drop = FALSE],
elev_im[list(x = stp0$df$x, y = stp0$df$y), drop = FALSE]
),
cbind(
land_use_f(dum1$x * 1000, dum1$y * 1000),
slope_im[list(x = dum1$x, y = dum1$y), drop = FALSE],
elev_im[list(x = dum1$x, y = dum1$y), drop = FALSE]
)
)
},
name = "covs0_spatial",
cache_dir = here("qcache")
)
## [1] "cached"
colnames(covs0_spatial) <- c("land", "slope", "elev")For the spatio-temporal case we need a more complicated lookup approach that takes into account the temporal dimension of the data. This is implemented in the function st_extract. First, we need to combine stp0 and dum1 into a unique object (so we don’t need to run st_extract twice which is extremely computationally intensive):
pts_combined <- rbind(stp0$df, dum1) |> setNames(c("x", "y", "time"))
pts_combined[, c("x", "y")] <- pts_combined[, c("x", "y")] * 1000 # Since we previously converted m to km
pts_combined <- st_as_sf(pts_combined, coords = c("x", "y"), crs = 3003)
pts_combined$time <- (as.Date("2023-01-01") + pts_combined$time) |> round_date(unit = "day")and reload the ERA-5 environmental data, taking their daily maxima value
env_var <- qcache(
{
tmp_dir <- tempdir()
unzip(here("data", "environmental-variables.zip"), exdir = tmp_dir)
nc_files <- list.files(
paste0(tmp_dir, "/environmental-variables"),
pattern = "\\.nc$",
full.names = TRUE
)
img_combined <- lapply(
nc_files,
function(x, var) {
img <- read_ncdf(x, var = var)
if (grepl("october.nc", x, fixed = TRUE)) {
img <- img[, , , 2, , drop = TRUE]
}
img
},
var = c("u10", "v10", "d2m", "t2m", "skt", "stl1", "stl2", "stl3", "stl4", "sp", "tp")
) |>
do.call(c, args = _)
unlink(nc_files); rm(nc_files); invisible(gc())
st_crs(img_combined) <- "OGC:CRS84"
st_warp(img_combined, crs = 3003)
},
name = "env_var",
cache_dir = here("qcache")
)
## [1] "cached"
env_var <- aggregate(env_var, max, by = "1 day")
invisible(gc())Now we can extract the ERA5-covariates values at these points
st_dimensions(env_var)[["time"]]$values <- seq.Date(
as.Date("2022-12-31"), by = "1 day", length.out = 365
)
env_cov_st <- st_extract(
x = env_var,
at = pts_combined,
time_column = "time"
)
env_cov_st <- st_drop_geometry(env_cov_st)[, -c(12)]and repeat the same operation for NDVI data
NDVI_cov_st <- st_extract(
x = NDVI_tidy,
at = pts_combined |> mutate(time = round_date(time, "month")),
time_column = "time",
)Let’s combine all together:
covs_complete <- cbind(
x = c(stp0$df$x, dum1$x),
y = c(stp0$df$y, dum1$y),
t = c(stp0$df$t, dum1$t),
covs0_spatial,
env_cov_st,
NDVI = NDVI_cov_st$NDVI
)
covs_complete$elev <- covs_complete$elev / 1000This is the object that will be passed to the modelling function which we can now define (together with several other auxiliary routines):
.counting.weights <- function(id, volumes) {
id <- as.integer(id)
fid <- factor(id, levels = seq_along(volumes))
counts <- table(fid)
w <- volumes[id] / counts[id]
w <- as.vector(w)
names(w) <- NULL
return(w)
}
.default.ncube <- function(X){
guess.ngrid <- floor((splancs::npts(X) / 2) ^ (1 / 3))
max(5, guess.ngrid)
}
.grid1.index <- function(x, xrange, nx) {
i <- ceiling(nx * (x - xrange[1]) / diff(xrange))
i <- pmax.int(1, i)
i <- pmin.int(i, nx)
i
}
.grid.index <- function(x, y, t, xrange, yrange, trange, nx, ny, nt) {
ix <- .grid1.index(x, xrange, nx)
iy <- .grid1.index(y, yrange, ny)
it <- .grid1.index(t, trange, nt)
return(list(ix = ix, iy = iy, it = it, index = as.integer((iy - 1) * nx + ix + (it - 1) * nx * ny)))
}
stppm <- function(X, formula, dummy_points, dati.interpolati,
ncube = NULL, obsvol, nxyt = NULL, local = FALSE,
verbose = TRUE) {
if (!inherits(X, c("stp", "stpm"))) {
stop("x should be either of class stp or stpm")
}
time1 <- Sys.time()
if (!is.null(ncube)) {
if (!is.numeric(ncube)) {
stop("ncube should be a numeric value")
} else {
if (ncube <= 0) {
stop("ncube should be ncube > 0")
}
}
}
# X è il processo osservato
X0 <- X
X <- X$df
nX <- nrow(X)
x <- X[, 1]
y <- X[, 2]
t <- X[, 3]
# definire dummy points come un df con x y e t
quad_p <- rbind(X[, 1:3], dummy_points)
xx <- quad_p[, 1]
xy <- quad_p[, 2]
xt <- quad_p[, 3]
win <- spatstat.geom::box3(
xrange = range(xx), yrange = range(xy),
zrange = range(xt)
)
if (is.null(ncube)) {
ncube <- .default.ncube(quad_p)
}
ncube <- rep.int(ncube, 3)
nx <- ncube[1]
ny <- ncube[2]
nt <- ncube[3]
if (is.null(nxyt)) nxyt <- nx * ny * nt
# cubevolume <- spatstat.geom::volume(win) / nxyt
cubevolume <- obsvol / nxyt
volumes <- rep.int(cubevolume, nxyt)
id <- .grid.index(
xx, xy, xt, win$xrange, win$yrange, win$zrange,
nx, ny, nt
)$index
w <- .counting.weights(id, volumes)
Wdat <- w[1:nX]
Wdum <- w[-c(1:nX)]
ndata <- nrow(X)
ndummy <- nrow(dummy_points)
# dati.interpolati sono i dati delle covariate
z <- c(rep(1, ndata), rep(0, nrow(dummy_points)))
y_resp <- z / w
dati.cov.marks <- data.frame(cbind(y_resp, w, quad_p, dati.interpolati))
suppressWarnings(mod_global <- try(gam(
as.formula(paste("y_resp",
paste(formula, collapse = " "),
sep = " "
)),
weights = w,
family = poisson,
data = dati.cov.marks
), silent = T))
summary(mod_global)
res_global <- coef(mod_global)
pred_global <- exp(predict(mod_global, newdata = dati.cov.marks[1:nX, ]))
if (local) {
nU <- dim(quad_p)[1]
h_x <- MASS::bandwidth.nrd(x)
h_y <- MASS::bandwidth.nrd(y)
h_t <- MASS::bandwidth.nrd(t)
localwt <- matrix(NA, nrow = nX, ncol = nU)
if (verbose) {
cat(paste(
"\n", "Computing Kernel Densities to the",
nX, "points", "\n", "\n"
))
}
for (j in 1:nX) {
if (verbose) {
spatstat.geom::progressreport(j, nX)
}
localwt[j, ] <- dnorm(xx - x[j], sd = h_x) *
dnorm(xy - y[j], sd = h_y) * dnorm(xt - t[j], sd = h_t)
}
a_s <- localwt * w
res_local <- matrix(NA, nrow = nX, ncol = length(mod_global$coefficients))
pred_local <- vector(length = nX)
if (verbose) {
cat(paste(
"\n", "Fitting local model to the", nX, "points",
"\n", "\n"
))
}
for (i in 1:nX) {
if (verbose) {
spatstat.geom::progressreport(i, nX)
}
suppressWarnings(mod_local <- try(glm(as.formula(paste("y_resp",
paste(formula, collapse = " "),
sep = " "
)), weights = a_s[i, ], family = poisson, data = dati.cov.marks), silent = T))
res_local[i, ] <- mod_local$coefficients
pred_local[i] <- exp(predict(mod_local, newdata = dati.cov.marks[i, ]))
}
res_local <- as.data.frame(res_local)
colnames(res_local) <- names(res_global)
}
time2 <- Sys.time()
if (local) {
list.obj <- list(
IntCoefs = res_global,
IntCoefs_local = res_local,
X = X0,
nX = ndata,
I = z,
y_resp = y_resp,
formula = formula,
l = as.vector(pred_global),
l_local = as.vector(pred_local),
mod_global = mod_global,
newdata = dati.cov.marks[1:ndata, ],
ncube = ncube,
time = paste0(round(as.numeric(difftime(
time1 = time2,
time2 = time1,
units = "sec"
)), 3), " sec")
)
class(list.obj) <- "locstppm"
} else {
list.obj <- list(
IntCoefs = res_global,
X = X0,
nX = ndata,
I = z,
y_resp = y_resp,
formula = formula,
l = as.vector(pred_global),
mod_global = mod_global,
newdata = dati.cov.marks[1:ndata, ],
ncube = ncube,
time = paste0(round(as.numeric(difftime(
time1 = time2,
time2 = time1,
units = "sec"
)), 3), " sec")
)
class(list.obj) <- "stppm"
}
# class(list.obj) <- "stppm"
return(list.obj)
}To conclude this step, we need to calculate the approximate volume of the 3D domain:
fires_sicily_ppp <- as.ppp(st_geometry(fires_sicily))
Window(fires_sicily_ppp) <- Window(land_use_tess)
volume0 <- area(
Window(fires_sicily_ppp |> spatstat.geom::rescale(1000))
) * abs(diff(range(stp0$df$t)))and the following is used to estimate the global model:
mod_global <- stppm(
X = stp0,
formula = ~
as.factor(land) + NDVI +
+ elev + slope
+ u10 + stl2
+ tp,
dummy_points = dum1,
dati.interpolati = covs_complete[3 + c(1, 2, 3, 4, 10, 13, 14, 15)],
obsvol = volume0, ncube = 5, nxyt = 87
)Check the parameters
summary(mod_global$mod_global)
##
## Family: poisson
## Link function: log
##
## Formula:
## y_resp ~ as.factor(land) + NDVI + +elev + slope + u10 + stl2 +
## tp
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -43.274352 0.864805 -50.039 < 2e-16 ***
## as.factor(land)2 -0.022468 0.048482 -0.463 0.643051
## as.factor(land)3 -0.205107 0.058849 -3.485 0.000492 ***
## as.factor(land)4 0.210173 0.182677 1.151 0.249931
## NDVI -4.790957 0.106466 -45.000 < 2e-16 ***
## elev 0.249990 0.032848 7.610 2.73e-14 ***
## slope 0.021782 0.001516 14.372 < 2e-16 ***
## u10 0.034866 0.006626 5.262 1.43e-07 ***
## stl2 0.127833 0.002840 45.005 < 2e-16 ***
## tp -54.544832 7.469982 -7.302 2.84e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## R-sq.(adj) = 0.179 Deviance explained = 38.1%
## UBRE = -0.61618 Scale est. = 1 n = 54781The following can be used to estimate the local model:
mod_local <- qcache(
{
stppm(
X = stp0,
formula = ~
as.factor(land) + NDVI +
+elev + slope
+ u10 + stl2
+ tp,
dummy_points = dum1,
dati.interpolati = covs_complete[3 + c(1, 2, 3, 4, 10, 13, 14, 15)],
obsvol = volume0, ncube = 5, nxyt = 87,
local = TRUE,
verbose = TRUE
)
},
name = "mod_local",
cache_dir = here("qcache")
)
## [1] "cached"Figure 13
The parameter estimates for the local model are summarised in the following plot:
plot_covs_3D <- function(
name,
bias,
main,
min = -Inf,
max = Inf,
val = NA,
mod = mod_local) {
id <- if (!is.na(val)) {
which(mod[["newdata"]][["land"]] == val)
} else {
seq_len(nrow(mod[["newdata"]]))
}
colvar <- mod[["IntCoefs_local"]][id, ][[name]]
id <- id[between(colvar, min, max)]
palette <- rev(colorRampPalette(
divergingx_hcl(
n = 11, l3 = 0, palette = "RdBu",
p3 = 0.8, p4 = 0.6
),
bias = bias
)(99))
x <- mod[["X"]][["df"]][["x"]][id]
y <- mod[["X"]][["df"]][["y"]][id]
t <- mod[["X"]][["df"]][["t"]][id]
scatter3D(
x = x,
y = y,
z = t,
theta = 60,
phi = 35,
col = palette,
ticktype = "detailed",
colvar = mod[["IntCoefs_local"]][id, name],
xlab = "Easting (Km)",
ylab = "Northing (Km)",
zlab = "Days from 1st January",
main = main
)
}if (!interactive()) {
par(mfrow = c(3, 3))
tmp <- mapply(
FUN = plot_covs_3D,
name = c(
"(Intercept)", "as.factor(land)2", "as.factor(land)3", "NDVI", "elev", "slope", "u10", "stl2", "tp"
),
bias = c(4, 0.85, 2.7, 5, 1.3, 0.9, 0.9, 0.2, 2.5),
main = c(
"Urban", "Agricultural Areas", "Forests", "NDVI", "Elevation", "Slope", "Wind Speed", "Temperature", "Precipitation"
),
min = c(-55, -1, -Inf, -5, -Inf, -Inf, -0.5, -Inf, -300),
max = c(-10, 1.25, 1, -3, 1, Inf, Inf, 0.25, 75),
val = c(1L, 2L, 3L, NA, NA, NA, NA, NA, NA)
)
rm(tmp)
}Figure 14
Global and local estimates of the intensity function:
if (!interactive()) {
# Intensity global
mark_int <- na.omit(mod_global$l)
ppp_gl <- ppp(
x = mod_global$X$df$x[-attr(mark_int, "na.action")],
y = mod_global$X$df$y[-attr(mark_int, "na.action")],
marks = mark_int,
window = mainland |> as.owin() |> spatstat.geom::rescale(1000)
)
int_global <- Smooth(ppp_gl, sigma = OS(unmark(ppp_gl)))
mark_int_local <- na.omit(mod_local$l_local)
ppp_l <- ppp(
x = mod_local$X$df$x[-attr(mark_int_local, "na.action")],
y = mod_local$X$df$y[-attr(mark_int_local, "na.action")],
marks = mark_int_local,
window = mainland |> as.owin() |> spatstat.geom::rescale(1000)
)
int_local <- Smooth(ppp_l, sigma = OS(unmark(ppp_l)))
int <- solist(int_global, int_local)
plot(
int,
equal.ribbon = TRUE,
main = "",
col = attr(spatstat.geom::colourmap(
grDevices::hcl.colors(100, "Viridis", rev = TRUE),
range = range(mark_int, mark_int_local)
),"stuff")$outputs,
main.panel = c("Global", "Local"),
panel.end = unmark(ppp_l),
panel.end.args = list(cex = 0.75)
)
}Figure 15
Comparison between the residuals in the global and local model.
# Focus on the mainland
id1 <- (!(fires_sicily_ppp$x < 1800000 | fires_sicily_ppp$y > 4263000))
# Remove a few points which displayed numerical problems in the estimate
id2 <- (
mod_local$IntCoefs_local$`tp` > c(-800)
| mod_local$IntCoefs_local$`stl2` < c(-0.2)
)
# Remove a few points with missing covariates due to numerical inaccuracies in the spatial extraction which implied NA in the estimates of lambda
id3 <- !is.na(mod_local$l_local)
id <- which(id1 & id2 & id3)
mod_local$newdata <- mod_local$newdata[id, ]
subset_ppp <- ppp(
x = mod_local$newdata$x,
y = mod_local$newdata$y,
marks = mod_local$newdata$t,
window = fires_sicily_ppp$window |> spatstat.geom::rescale(1000)
)
# Spatiotemporal kernel intensity
st_kernel <- spattemp.density(subset_ppp, tres = 128)
# Extract the intensity values. The warnings you might see are simply due to
# numerical rounding problems at the boundary of the time interval.
cc <- spattemp.slice(st_kernel, tt = mod_local$newdata$t)$z
cc_dens = mapply(
FUN = function(x, y, cc) {
cc[list(x = x, y = y)]
},
x = mod_local$newdata$x,
y = mod_local$newdata$y,
cc = cc
)
# We need to remove such points at the boundary of the time interval
id_no <- vapply(cc_dens, identical, logical(1), y = numeric(0)) |> which()
nX <- length(unlist(cc_dens))
dt <- data.frame(
pred = c(
mod_local$l[id][-id_no], mod_local$l_local[id][-id_no]
),
smo = c(unlist(cc_dens) * nX, unlist(cc_dens) * nX),
type = rep(c("Global", "Local"), each = length(unlist(cc_dens)))
)if (!interactive()) {
ggplot(dt, aes(x = pred, y = smo)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", linewidth = 1, col = 2) +
facet_wrap(~type, scales = "free") +
scale_x_continuous(transform = "log10") +
scale_y_continuous(transform = "log10") +
theme_light() +
theme(
strip.text = element_text(size = 11, colour = "black", face = "bold"),
axis.text = element_text(size = 8.5),
axis.title = element_text(size = 11.5)
) +
labs(x = "Parametric intensity", y = "Kernel intensity")
}Figure 16
Density distribution of the residuals
residui_global <- unlist(cc_dens) * nX - mod_local$l[id][-id_no]
residui_local <- unlist(cc_dens) * nX - mod_local$l_local[id][-id_no]
residui <- data.frame(
residuals = c(residui_global, residui_local),
model = rep(c("Global", "Local"), each = length(residui_local))
)if (!interactive()) {
ggplot(data=residui, aes(x=residuals, group=model, fill=model)) +
geom_density(adjust=1.5, alpha=.4) +
theme_minimal() +
theme(
axis.text = element_text(size = 10)
) +
labs(x = "Residual", y = "Density", fill = "Model Type")
}Figure 17
Spatio-temporal residuals of the local model:
color_palette_on_zero <- rev(colorRampPalette(divergingx_hcl(
n = 11, l3 = 0,
palette = "RdBu",
p3 = 0.8, p4 = 0.6
), bias = 1.5)(99))
if (!interactive()) {
scatter3D(
x = mod_local$newdata$x[-id_no],
y = mod_local$newdata$y[-id_no],
z = mod_local$newdata$t[-id_no],
col = color_palette_on_zero,
theta = 60, phi = 35,
colvar = residui_local,
ticktype = "detailed",
xlab = "Northing (Km)", ylab = "Easting (Km)", zlab = "Days from January 1st"
)
}